using chc 6
present CV LDAs, have testing/training if need be
graphed LDAs are all the points.

MANOVA for species level

fit_full_species_man$pca_summary 
## Importance of first k=7 (out of 35) components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     1.0413 0.7144 0.51001 0.38099 0.34918 0.30663 0.28705
## Proportion of Variance 0.3871 0.1822 0.09287 0.05183 0.04354 0.03357 0.02942
## Cumulative Proportion  0.3871 0.5694 0.66225 0.71408 0.75762 0.79119 0.82061
summary(manova(as.matrix(data[,4:cols]) ~ hostRace * sex *site, data = data), test = "Wilks")
##                    Df   Wilks approx F num Df den Df    Pr(>F)    
## hostRace            4 0.04507   53.399     28 1097.5 < 2.2e-16 ***
## sex                 1 0.71370   17.421      7  304.0 < 2.2e-16 ***
## site                8 0.26950    8.085     56 1642.4 < 2.2e-16 ***
## hostRace:sex        4 0.80887    2.375     28 1097.5 8.147e-05 ***
## hostRace:site       1 0.96866    1.405      7  304.0  0.202640    
## sex:site            7 0.77071    1.663     49 1547.8  0.002979 ** 
## hostRace:sex:site   1 0.97168    1.266      7  304.0  0.266990    
## Residuals         310                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PCA for Species Level

LDA for Species level

##            Reference
## Prediction  Cingulata Cornivora Mendax pom Zepheria
##   Cingulata        29         0      0   0        0
##   Cornivora         1         5      0   1        0
##   Mendax            0         0     60   8        2
##   pom               0         0      1 214        0
##   Zepheria          0         0      0   1       15
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##   9.584570e-01   9.213168e-01   9.312796e-01   9.771049e-01   6.646884e-01 
## AccuracyPValue  McnemarPValue 
##   2.679678e-40            NaN

lda for species + sex + site

##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##   5.164179e-01   4.800429e-01   4.614714e-01   5.710717e-01   1.432836e-01 
## AccuracyPValue  McnemarPValue 
##   4.176925e-58            NaN

Manova for pomonella

summary(manova(as.matrix(data[,4:cols]) ~ hostRace * site * sex, data = data), test = "Wilks")
##                    Df   Wilks approx F num Df den Df    Pr(>F)    
## hostRace            1 0.70245  12.2841      7 203.00 4.504e-13 ***
## site                4 0.37641   8.1524     28 733.35 < 2.2e-16 ***
## sex                 1 0.48494  30.8012      7 203.00 < 2.2e-16 ***
## hostRace:site       2 0.72182   5.1338     14 406.00 6.340e-09 ***
## hostRace:sex        1 0.95257   1.4440      7 203.00  0.189418    
## site:sex            4 0.77940   1.8745     28 733.35  0.004284 ** 
## hostRace:site:sex   1 0.92648   2.3014      7 203.00  0.028089 *  
## Residuals         209                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

LDA for Pomonella

##           Reference
## Prediction Apple Haw
##      Apple    80  22
##      Haw      26  96
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##   7.857143e-01   5.693688e-01   7.261263e-01   8.375728e-01   5.267857e-01 
## AccuracyPValue  McnemarPValue 
##   8.355647e-16   6.650055e-01

LDA predicting host and sex for pomonella

##               Reference
## Prediction     Apple_Female Apple_Male Haw_Female Haw_Male
##   Apple_Female           51          5         17        2
##   Apple_Male              1         21          1        8
##   Haw_Female             19          1         44        5
##   Haw_Male                0          8          4       37
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##   6.830357e-01   5.662230e-01   6.177391e-01   7.434095e-01   3.169643e-01 
## AccuracyPValue  McnemarPValue 
##   3.310270e-29   5.581411e-01

MANOVA for Zepheria

# no interaction because missing males at MtPleasant
summary(manova(as.matrix(data[,4:cols]) ~ sex + site, data = data), test = "Wilks")
##           Df   Wilks approx F num Df den Df Pr(>F)
## sex        1 0.74708  0.67708      5     10 0.6507
## site       1 0.83749  0.38808      5     10 0.8461
## Residuals 14

LDA for Zepheria

##                  Reference
## Prediction        Zepheria_Female Zepheria_Male
##   Zepheria_Female               5             3
##   Zepheria_Male                 4             5
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##      0.5882353      0.1793103      0.3292472      0.8155630      0.5294118 
## AccuracyPValue  McnemarPValue 
##      0.4062810      1.0000000

LDA for sex +site zepheria

##                              Reference
## Prediction                    Zepheria_Female_EastLansing
##   Zepheria_Female_EastLansing                           2
##   Zepheria_Female_MtPleasant                            0
##   Zepheria_Male_EastLansing                             4
##                              Reference
## Prediction                    Zepheria_Female_MtPleasant
##   Zepheria_Female_EastLansing                          2
##   Zepheria_Female_MtPleasant                           1
##   Zepheria_Male_EastLansing                            0
##                              Reference
## Prediction                    Zepheria_Male_EastLansing
##   Zepheria_Female_EastLansing                         1
##   Zepheria_Female_MtPleasant                          5
##   Zepheria_Male_EastLansing                           2
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##     0.29411765    -0.05699482     0.10313551     0.55958272     0.47058824 
## AccuracyPValue  McnemarPValue 
##     0.95792652     0.03207164

MANOVA for Mendax

summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")
##           Df   Wilks approx F num Df den Df    Pr(>F)    
## sex        1 0.57727   6.1024      6     50 7.541e-05 ***
## site       2 0.10716  17.1232     12    100 < 2.2e-16 ***
## sex:site   2 0.74705   1.3082     12    100    0.2258    
## Residuals 55                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

LDA for Mendax

##                Reference
## Prediction      Mendax_Female Mendax_Male
##   Mendax_Female            20          15
##   Mendax_Male              13          13
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##     0.54098361     0.07072905     0.40849889     0.66935590     0.54098361 
## AccuracyPValue  McnemarPValue 
##     0.55241652     0.85010674

LDA for sex +site mendax

##                         Reference
## Prediction               Mendax_Female_Fenville Mendax_Female_OtisLake
##   Mendax_Female_Fenville                      7                      1
##   Mendax_Female_OtisLake                      1                      5
##   Mendax_Female_Sewanee                       0                      2
##   Mendax_Male_Fenville                        6                      0
##   Mendax_Male_OtisLake                        0                      1
##   Mendax_Male_Sewanee                         0                      0
##                         Reference
## Prediction               Mendax_Female_Sewanee Mendax_Male_Fenville
##   Mendax_Female_Fenville                     0                    5
##   Mendax_Female_OtisLake                     1                    0
##   Mendax_Female_Sewanee                      4                    1
##   Mendax_Male_Fenville                       0                    5
##   Mendax_Male_OtisLake                       1                    1
##   Mendax_Male_Sewanee                        4                    0
##                         Reference
## Prediction               Mendax_Male_OtisLake Mendax_Male_Sewanee
##   Mendax_Female_Fenville                    0                   0
##   Mendax_Female_OtisLake                    1                   2
##   Mendax_Female_Sewanee                     2                   3
##   Mendax_Male_Fenville                      0                   1
##   Mendax_Male_OtisLake                      1                   1
##   Mendax_Male_Sewanee                       2                   3
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##    0.409836066    0.283523654    0.285504382    0.543223627    0.229508197 
## AccuracyPValue  McnemarPValue 
##    0.001280557            NaN

MANOVA for cingulata

summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")
##           Df   Wilks approx F num Df den Df   Pr(>F)   
## sex        1 0.80383   1.2813      4     21 0.308955   
## site       2 0.39036   3.1529      8     42 0.006927 **
## sex:site   2 0.66131   1.2059      8     42 0.319040   
## Residuals 24                                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

LDA for Cingulata

##                   Reference
## Prediction         Cingulata_Female Cingulata_Male
##   Cingulata_Female                8              6
##   Cingulata_Male                  7              9
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##      0.5666667      0.1333333      0.3742735      0.7453925      0.5000000 
## AccuracyPValue  McnemarPValue 
##      0.2923324      1.0000000

LDA for sex +site cingulata

##                             Reference
## Prediction                   Cingulata_Female_Fenville
##   Cingulata_Female_Fenville                          1
##   Cingulata_Female_SouthBend                         0
##   Cingulata_Female_Urbana                            0
##   Cingulata_Male_Fenville                            1
##   Cingulata_Male_SouthBend                           0
##   Cingulata_Male_Urbana                              1
##                             Reference
## Prediction                   Cingulata_Female_SouthBend Cingulata_Female_Urbana
##   Cingulata_Female_Fenville                           3                       0
##   Cingulata_Female_SouthBend                          2                       0
##   Cingulata_Female_Urbana                             4                       1
##   Cingulata_Male_Fenville                             0                       0
##   Cingulata_Male_SouthBend                            0                       0
##   Cingulata_Male_Urbana                               1                       1
##                             Reference
## Prediction                   Cingulata_Male_Fenville Cingulata_Male_SouthBend
##   Cingulata_Female_Fenville                        0                        4
##   Cingulata_Female_SouthBend                       0                        1
##   Cingulata_Female_Urbana                          0                        1
##   Cingulata_Male_Fenville                          0                        1
##   Cingulata_Male_SouthBend                         0                        1
##   Cingulata_Male_Urbana                            0                        2
##                             Reference
## Prediction                   Cingulata_Male_Urbana
##   Cingulata_Female_Fenville                      0
##   Cingulata_Female_SouthBend                     2
##   Cingulata_Female_Urbana                        0
##   Cingulata_Male_Fenville                        2
##   Cingulata_Male_SouthBend                       0
##   Cingulata_Male_Urbana                          0
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##     0.17241379     0.04000000     0.05845608     0.35774755     0.34482759 
## AccuracyPValue  McnemarPValue 
##     0.98826571            NaN

MANOVA for Cornivora

There is only one site

# no site because only sampled at one site.
summary(manova(as.matrix(data[,4:cols]) ~ sex, data = data), test = "Wilks")
##           Df   Wilks approx F num Df den Df Pr(>F)
## sex        1 0.54944  0.82004      2      2 0.5494
## Residuals  3

LDA for cornivora

##                   Reference
## Prediction         Cornivora_Female Cornivora_Male
##   Cornivora_Female                0              2
##   Cornivora_Male                  0              2
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##     0.50000000     0.00000000     0.06758599     0.93241401     1.00000000 
## AccuracyPValue  McnemarPValue 
##     1.00000000     0.47950012

lines of diversification

euclidian distance

##                  Apple_Female Apple_Male Cingulata_Female Cingulata_Male
## Apple_Female        0.0000000                                           
## Apple_Male          1.1052241  0.0000000                                
## Cingulata_Female    1.4738114  1.6036342        0.0000000               
## Cingulata_Male      1.4051571  1.6128575        0.1130547      0.0000000
## Cornivora_Female    1.8164505  1.6514454        0.4825512      0.5956054
## Cornivora_Male      2.3568964  1.7473686        1.3312107      1.4407462
## Haw_Female          0.5292825  1.5234184        1.9508336      1.8698619
## Haw_Male            0.6482415  0.5290546        1.6459586      1.6184269
## Mendax_Female       1.1742028  2.0634150        1.2374012      1.1243501
## Mendax_Male         1.0269548  1.8389388        1.0200228      0.9078379
## Zepheria_Female     2.2810030  3.3164335        2.4854319      2.3757562
## Zepheria_Male       2.3942415  3.4468577        2.6586999      2.5484713
##                  Cornivora_Female Cornivora_Male Haw_Female  Haw_Male
## Apple_Female                                                         
## Apple_Male                                                           
## Cingulata_Female                                                     
## Cingulata_Male                                                       
## Cornivora_Female        0.0000000                                    
## Cornivora_Male          0.8731903      0.0000000                     
## Haw_Female              2.3257992      2.8860880  0.0000000          
## Haw_Male                1.8371458      2.1364704  1.0022950 0.0000000
## Mendax_Female           1.7199319      2.5456134  1.3364646 1.7570533
## Mendax_Male             1.5002006      2.3102051  1.2814072 1.5661386
## Zepheria_Female         2.9567427      3.8164525  2.1620829 2.9257376
## Zepheria_Male           3.1317505      3.9891585  2.2408291 3.0416391
##                  Mendax_Female Mendax_Male Zepheria_Female Zepheria_Male
## Apple_Female                                                            
## Apple_Male                                                              
## Cingulata_Female                                                        
## Cingulata_Male                                                          
## Cornivora_Female                                                        
## Cornivora_Male                                                          
## Haw_Female                                                              
## Haw_Male                                                                
## Mendax_Female        0.0000000                                          
## Mendax_Male          0.2497539   0.0000000                              
## Zepheria_Female      1.3207469   1.5704754       0.0000000              
## Zepheria_Male        1.4794064   1.7289277       0.1809107     0.0000000

mds plot

centroids

centroid slopes

** should I calc slope for each host race pair? or only two groups as shown below?

pomonella == y ~ -0.8633 + 1.4937 * LD2.mean
rest == y ~ 1.8431 + 1.5561 * LD2.mean

these slopes are very similar. Just different intercepts. relationship btwn axes is the same for all species, but the pomonella host races have greater LD2 values compared to the other species. Also female to male differentiation flows along the same slope except the relationship is reversed in zepheria and cingulata. Cingulata, zepheria, and mendax have small sex diff. Pomonella has great sex differences. Cornivora look to have great sex differences but further sampling will be needed to confrim.

## 
## Call:
## lm(formula = LD1.mean ~ LD2.mean, data = .)
## 
## Residuals:
##        1        2        3        4 
##  0.40664 -0.10011 -0.08027 -0.22626 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -0.8633     0.2246  -3.843   0.0615 .
## LD2.mean      1.4937     0.5724   2.609   0.1208  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3413 on 2 degrees of freedom
## Multiple R-squared:  0.7729, Adjusted R-squared:  0.6594 
## F-statistic: 6.808 on 1 and 2 DF,  p-value: 0.1208
## 
## Call:
## lm(formula = LD1.mean ~ LD2.mean, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.48301 -0.15745  0.06621  0.23799  0.31705 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.8431     0.1729  10.659 4.02e-05 ***
## LD2.mean      1.5561     0.1606   9.686 6.94e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3161 on 6 degrees of freedom
## Multiple R-squared:  0.9399, Adjusted R-squared:  0.9299 
## F-statistic: 93.83 on 1 and 6 DF,  p-value: 6.944e-05

mcmc sampling ???

calc angle

theta = tan_inverse (slope)

slope = (y1 - y0)/(x1 - x0)
theta = atan(slope)
I think the slope is given below in radians? y1, x1= always the female and y0 an x0 are always males.it should’t matter though…

# apple
ang = centroids %>% filter(host_Species == "Apple")
ang
##   host_Species    sex  LD1.mean   LD2.mean
## 1        Apple Female -0.115682  0.2282819
## 2        Apple   Male -1.160544 -0.1319735
y1 = ang$LD2.mean[1] # female y
y0 = ang$LD2.mean[2] #male y
x1 = ang$LD1.mean[1] # female x
x0 = ang$LD1.mean[2] #male x
slope = (y1 - y0)/(x1 - x0)
slope
## [1] 0.3447876
theta = atan(slope)
theta
## [1] 0.3320237
apple = cbind(slope,theta)
apple = as.data.frame(apple)
apple$host = "apple"
apple
##       slope     theta  host
## 1 0.3447876 0.3320237 apple
# haw
ang = centroids %>% filter(host_Species == "Haw")
ang
##   host_Species    sex   LD1.mean  LD2.mean
## 1          Haw Female  0.1114762 0.7063395
## 2          Haw   Male -0.7638430 0.2180654
y1 = ang$LD2.mean[1] # female y 
y0 = ang$LD2.mean[2] #male y 
x1 = ang$LD1.mean[1] # female x 
x0 = ang$LD1.mean[2] #male x
slope = (y1 - y0)/(x1 - x0)
theta = atan(slope)

haw = cbind(slope,theta)
haw = as.data.frame(haw)
haw$host = "haw"
haw
##       slope     theta host
## 1 0.5578241 0.5088304  haw
# zepheria
ang = centroids %>% filter(host_Species == "Zepheria")
ang
##   host_Species    sex LD1.mean     LD2.mean
## 1     Zepheria Female 2.153432 -0.004308698
## 2     Zepheria   Male 2.276472  0.128317729
y1 = ang$LD2.mean[1] # female y 
y0 = ang$LD2.mean[2] #male y 
x1 = ang$LD1.mean[1] # female x 
x0 = ang$LD1.mean[2] #male x 
slope = (y1 - y0)/(x1 - x0)
theta = atan(slope)

Zepheria = cbind(slope,theta)
Zepheria = as.data.frame(Zepheria)
Zepheria$host = "Zepheria"
Zepheria
##      slope     theta     host
## 1 1.077911 0.8228755 Zepheria
# mendax
ang = centroids %>% filter(host_Species == "Mendax")
ang
##   host_Species    sex  LD1.mean   LD2.mean
## 1       Mendax Female 0.8877174 -0.3815885
## 2       Mendax   Male 0.6495075 -0.4566423
y1 = ang$LD2.mean[1] # female y 
y0 = ang$LD2.mean[2] #male y
x1 = ang$LD1.mean[1] # female x 
x0 = ang$LD1.mean[2] #male x 
slope = (y1 - y0)/(x1 - x0)
theta = atan(slope)

Mendax = cbind(slope,theta)
Mendax = as.data.frame(Mendax)
Mendax$host = "Mendax"
Mendax
##       slope     theta   host
## 1 0.3150745 0.3052286 Mendax
# cingulata
ang = centroids %>% filter(host_Species == "Cingulata")
ang
##   host_Species    sex     LD1.mean  LD2.mean
## 1    Cingulata Female -0.002386062 -1.241168
## 2    Cingulata   Male  0.079528433 -1.163249
y1 = ang$LD2.mean[1] # female y 
y0 = ang$LD2.mean[2] #male y
x1 = ang$LD1.mean[1] # female x 
x0 = ang$LD1.mean[2] #male x 
slope = (y1 - y0)/(x1 - x0)
theta = atan(slope)

Cingulata = cbind(slope,theta)
Cingulata = as.data.frame(Cingulata)
Cingulata$host = "Cingulata"
Cingulata
##       slope     theta      host
## 1 0.9512238 0.7604056 Cingulata
# Cornivora
ang = centroids %>% filter(host_Species == "Cornivora")
ang
##   host_Species    sex  LD1.mean  LD2.mean
## 1    Cornivora Female -0.353129 -1.572582
## 2    Cornivora   Male -1.170673 -1.879313
y1 = ang$LD2.mean[1] # female y 
y0 = ang$LD2.mean[2] #male y
x1 = ang$LD1.mean[1] # female x 
x0 = ang$LD1.mean[2] #male x 
slope = (y1 - y0)/(x1 - x0)
theta = atan(slope)

Cornivora = cbind(slope,theta)
Cornivora = as.data.frame(Cornivora)
Cornivora$host = "Cornivora"
Cornivora
##       slope     theta      host
## 1 0.3751856 0.3589334 Cornivora
rbind(apple, haw, Zepheria, Mendax, Cingulata, Cornivora)
##       slope     theta      host
## 1 0.3447876 0.3320237     apple
## 2 0.5578241 0.5088304       haw
## 3 1.0779109 0.8228755  Zepheria
## 4 0.3150745 0.3052286    Mendax
## 5 0.9512238 0.7604056 Cingulata
## 6 0.3751856 0.3589334 Cornivora

site centroids